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Mh Abstract: We analyze a quantum reservoir engineering method, originally introduced by Sar- 

lette et al [Phys. Rev. Lett. 107, 010402 (2011)], for the stabilization of non-classical field states 
in high quality cavities. We generalize the method to the protection of mesoscopic entangled field 

F^ states shared by two non-degenerate field modes. The reservoir is made up of a stream of atoms 

O^ undergoing successive composite interactions with the cavity, each combining resonant with non- 

'""' resonant parts. We get a detailed insight into the competition between the engineered reservoir 

^_^ and decoherence. We show that the operation is quite insensitive to experimental imperfections 

^ and that it could thus be implemented in the near future, either in the context of microwave 

On Cavity Quantum Electrodynamics or in that of circuit-QED. 
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en 1 Introduction 

O 

CN Nonclassical electromagnetic field states are extremely important, both for a fundamental under- 

. . standing of the quantum properties of light and for their possible use in practical applications. 

^ For instance, squeezed states (SS) have fluctuations of one of their quadratures below those of 

the vacuum state, or of a classical coherent state [1^. They lead thus to interesting methods for 
high-precision measurements and metrology [2]. They are for instance planned to be used for 
C^ reducing the noise of the gravitational wave interferometers below the standard quantum limit [3^ . 

Mesoscopic field state superpositions (MESS) are also the focus of an intense interest. They 
involve a quantum superposition of two quasi-classical coherent components with different com- 
plex amplitudes. These counter-intuitive states bridge the gap between the quantum and the 
classical worlds and shed light onto the decoherence process responsible for the conspicuous lack 
of superpositions at our scale [4] . 

Einally, entangled superpositions of mesoscopic states (ESMS) shared by several field modes are 
even more intriguing. They violate generalized Bell inequalities [5 , illustrating the fundamentally 
non-local nature of quantum phyics. However, their non-local character is rapidly erased by a 
fast decoherence process [6 , driving them back into a statistical mixture that can be undestood 
in terms of a classical local hidden variable model. This interplay of decoherence and nonlocality 
opens fascinating perspectives for exploring the limits of the quantum. 
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Figure 1: (a-d) Wigner functions of nonclassical field states e~*^^^^|a) generated by propagation 
of an initial coherent state through a Kerr medium, (a) 2-component MFSS given by Eq. ([2| for 
^KlK = 7r/2; (b) 3-component MFSS for tx^K = 7r/3; (c) 'banana'-state, for tx^K = 0.28; and 
(d) squeezed state, for tKjK = 0.08 <C tt. (e-h): similar states stabilized, despite decoherence, 
by the atomic reservoir as explained in the remainder of the paper. Frame (e) corresponds to the 
reference two-component MFSS most lengthily discussed in the rest of the paper. 



In principle, the SS and MFSS could be simply prepared in the optical domain by letting 
a coherent laser pulse propagate in a non-linear medium, whose index of refraction is a linear 
function of the light pulse intensity (Kerr medium) [T. The field evolves from initial coherent 
state \a) under the action of the Kerr Hamiltonian H^^ 



H 



K 



Ck N + 7K N^ 



(1) 



Here N is the photon number operator, (k is proportional to the linear index and 7^ is the 
Kerr frequency describing the strength of the non-linearity. In the following, we use units such 
that h = 1. Note that the collisional interaction Hamiltonian for an N-atom sample in a tightly 
confining potential or in an optical lattice is similar to H^ [8^. 

Depending on the interaction time tx, the final state e~*^^^^|a) spans a number of nonclas- 
sical forms m Section 7.2], including: 

(i) squeezed states for txjK <C tt; 

(ii) states with 'banana'-shaped Wigner function for slightly larger tx^K] 

(iii) mesoscopic field state superpositions \ka) with k equally spaced components for tx^K = 
7i/k [9]. 



(iv) in particular, a MFSS of two coherent states with opposite amplitudes: 

\ca) = {\a)^i\-a))/V2, 



(2) 



with a = ae *'^^^^, for txjK 



2 • 



The top panels on Figure IT] present the Wigner functions of the states (i)-(iv) for a mean photon 
number laP = 2.7. 



This preparation method by a deterministic unitary evolution is simple in its principle, but its 
implementation is extraordinarily difficult for propagating light fields, due to the weakness of the 
Kerr nonlinearity [10 . 

Other methods for the production of these nonclassical states have been proposed or realized 
in the context of trapped ions [TT1|12] or Cavity Quantum Electrodynamics (CQED) [4) [T3 l [M ) [T5 l 
[l6l[T7l[^, 19 . Both systems implement the 'spin-spring' model, the simplest nontrivial quantum 
situation of a two-level system coupled to a harmonic oscillator, embodied by the harmonic motion 
of the ion or by a single field mode. The proposed nonclassical state production methods are either 
deterministic or rely on a detection- conditioned scheme. The latter expand the possibilities of the 
former by applying a measurement operation after a unitary evolution towards an intermediate 
target state. Measurement back-action generates different final states conditioned by the stochastic 
detection outcome [20]. In the microwave CQED context, detection-conditioned preparation of 
MESS and ESMS can be achieved by the dispersive interaction of an initial coherent field state 
with a two-level atom, initially prepared in a state superposition, followed by the detection of the 
atomic state in an appropriate basis [4ll2T]. 

All these preparation techniques do not solve, however, the problem of stabilizing ("protecting") 
a selected nonclassical state for long times in spite of the unavoidable coupling of the system S 
to its environment £. Reservoir engineering can be used to stabilize target quantum states by 
strongly coupling S to an "engineered" environment, or reservoir 7^, a large quantum system 
with many degrees of freedom. The reservoir is designed so that, when acting alone, it drives <S, 
whatever its intial state, towards a unique target 'pointer state', a stable state of S coupled to 
71 [22l [23]. The state of S remains close to this pointer state even in the presence of £^ provided 
IZ is more strongly coupled to S than S. An engineered reservoir thus achieves much more than 
the preparation of a target state. It effectively stabilizes the system close to it for arbitrarily long 
times. 

Reservoir engineering is experimentally challenging. Reservoirs made up of lasers and mag- 
netic fields for trapped- ion oscillators have been proposed [24l [191 ES] ^^^ demonstrated [26 . 
Recently, a reservoir has been used to generate entanglement of spin states of macroscopic atomic 
ensembles [27] . 

In the context of CQED, the reservoir may be a stream of atoms interacting with the trapped 
field. An early proposal [28 relied on the so-called 'trapping state conditions' for the micro- 
maser [29] , which require a very fine tuning of the parameters and can only be properly achieved 
in the case of a zero-temperature environment. Reservoirs composed of atoms in combination with 
external fields have also been proposed to stabilize one-mode squeezed states [30] and two-mode 
squeezed vacuum entanglement [31 . 

In [32 , we proposed a robust reservoir engineering method for CQED. It generates and sta- 
bilizes nonclassical states of a single mode of the radiation field, including SS and MFSS. The 
reservoir is made up of a stream of 2-level atoms, each prepared in a coherent superposition of its 
basis states. They interact one at a time with the field according to the Jaynes-Cummings model 
before being discarded, a procedure reminiscent of the "reset" operation performed in other con- 
texts [33l[26]. The key idea is to use a tailored composite interaction of each atom with the field: 
dispersive, then resonant, then dispersive again. The pointer states of this composite interaction 
are precisely those, e~'^^^^^\a), resulting from the action of a Kerr Hamiltonian acting upon an 
initially coherent state. 

This method is quite general and could be implemented in a variety of CQED settings, par- 
ticularly in the active context of circuit QED [34 or in that of microwave CQED, with circular 
Rydberg atoms and superconducting Fabry Perot cavities. For the sake of definiteness, we shall 
focus in this paper on the microwave CQED case, and particularly on the current ENS CQED 
experiment whose scheme is depicted on Fig. [2] The bottom panels of Fig. [l] present the results 
of numerical simulations of the ENS experiment, with interaction parameters chosen to reproduce 
the states generated by the Kerr Hamiltonian (top panels). 

The present paper is intended to provide an in-depth description of this single-mode reservoir 
engineering procedure, with a detailed analysis of the physical mechanism of state stabilization. We 
discuss also the competition between the engineered reservoir and the ordinary cavity enviroment. 



giving simple insights into the finite final fidelity of the prepared state. 

We finally extend the scheme proposed in [32^ to the stabilization of entangled superpositions 
of mesoscopic states of two field modes. The atoms of the reservoir undergo a tailored interaction 
with two modes of the same cavity, combining dispersive and resonant parts for each mode. This 
proposal opens interesting perspectives for studying the interplay between entanglement, non- 
locality and decoherence in the context of mesoscopic quantum states. 

The paper is organized as follows. We consider the single- mode case for most of the paper 
and extend it to two modes in the last Section. Section [2] describes the experimental scheme and 
the principle of the method. Section [3] discusses, as a building block for the next Sections, how 
a stream of atoms resonant with one field mode stabilizes approximately a coherent field state. 
Section [4] introduces the composite interaction: non-resonant, resonant and non-resonant again. 
In this Section, we treat the non-resonant interactions in the dispersive regime of a large atom- 
cavity detuning. We thus get a simple qualitative insight into the mechanism generating non- 
classical states. Section [5] details the more realistic case of intermediate atom-cavity detuning. 
We show that the main features of Section [4] are recovered, exhibiting the robustness of the 
method. Section [6] analyzes the effect of decoherence due to the cavity damping and imprecisions 
on the experimental parameters. We find that the method is also robust against realistically large 
imperfections. Section [7| finally presents the stabilization of two- mode ESMS. 

2 General description 




Figure 2: Scheme of the current ENS CQED experiment. 

The scheme of the ENS experiment is depicted on Fig. |2] (see [4l[T5] for details). A microwave 
field of frequency ujc is trapped in the superconducting cavity C (damping time Tc = 65 ms). 
Atoms are sent one after the other through C. They cross its standing wave gaussian mode at a 
constant, adjustable velocity v. The mode interacts with the transition between the two atomic 
levels l^f) and |e) (circular Rydberg states with principal quantum numbers 50 and 51). A static 
electric field applied across the cavity mirrors is used to adjust the atomic transition frequency cjq 
and hence the atom-cavity detuning 5 = ujq — ujc ^ ^c via the Stark effect. Varying the electric 
field during the atom-field interaction makes it possible to engineer the detuning profile S{t). Zero 
and small S values are used for the resonant and non-resonant parts of the interaction. Making S 
very large allows us to effectively turn off the atom-field interaction. 

We describe the atom and field states in a frame rotating at frequency ujc. The atoms are 
prepared in state |^) in 5, by a time-resolved laser and radiofrequency excitation of a velocity- 
selected thermal rubidium atomic beam. Before entering C, the atoms are prepared in a coherent 
superposition of \g) and |e) in the low-quality cavity i?i, driven by a classical microwave source 
("first Ramsey zone") at frequency cjc- Without loss of generality, we can choose the phase refer- 
ence for all atoms so that they enter the cavity in the initial state \ug,t) = cos{u/2)\g) +sin('u/2)|e) 
with u > 0. In a Bloch sphere representation with |e) at the north pole, \ug,t) corresponds to a 
vector at an angle u with the north-south vertical axis. 

A second classical microwave pulse in the second Ramsey zone R2 is followed by a detection in 
the {|e), 1^)} basis in the field- ionization detector D. This amounts to a projective measurement of 
the atomic state at the exit of C, in a basis that can be chosen arbitrarily by properly setting the 



microwave pulse in R2. For the engineered reservoir operation, the result of this final atomic state 
detection is irrelevant. Detection results are however necessary in other phases of the experiment. 
In particular, they will be used to reconstruct the field state generated by the engineered reservoir, 
using a method described in [15 . 

Let us first consider atom-cavity interaction for an atom that crosses cavity axis at t = 0. 
The atom-field interaction is ruled by the Jaynes-Cummings Hamiltonian Hjc- Neglecting far 
off-resonant terms (rotating wave approximation, negligible interaction with other cavity modes), 
we have: 

Hjc = ^(|e)(e| - |5)(5|) + i^(|5)(e|at - \e){g\e.) , (3) 

where Q{s) is the atom-cavity coupling strength (vacuum Rabi frequency) at position s = vt along 
the atomic trajectory; a is the photon annihilation operator in C. The photon number operator 
N = a'f'a = ^^ n \n){n\ defines the Fock states basis {\n)}. 

The coupling strength Q{s) is determined by the atomic transtion parameters and by the cavity 
mode geometry. It writes here ft{s) = Qq e~^ 1^ , with Vl^l2i\: = 50 kHz and w = 6 mm for the 
ENS setup. To get a finite total interaction duration T, we assume that the coupling cancels when 
1^1 > 1.5 w. The total interaction time of the atom with the field is thus T = 3w/v. We have 
checked numerically that this truncation of interaction time has negligible effects. 

The evolution operator, or propagator U associated to Hjc expresses the transformation that 
the joint atom-cavity state undergoes during interaction. The Schrodinger equation for U, starting 
at the initial time t = to is: 



V{t) = -i Hjcit) V{t) with U(to) = I , (4) 



where I is the identity operator. We note Ut the propagator obtained by integration of Eq. Q 
over one full atom-cavity interaction, that lasts from — T/2 to T/2. 

We represent the action of Ut over the field state by the operators M^^ and M^^, such that: 

UT(kat>|^» = |5)My-|^) + |e)My-|^), 

for any pure initial field state |?/^). Tracing over the final atomic state, the modification of the field 
density operator due to the interaction with the atom is thus finally given by the Kraus map [35j 

p -^ M^^pM^^ ^ + M^^pM^^ ^ . (5) 

For the reservoir action, we let a stream of atoms consecutively interact with the field and 
always use the same parameter set (detuning profile, atom velocity v and initial state \ug,t) = 
cos{u/2)\g) -\- sm{u/2)\e)). Thus each atom affects the field according to Eq. (Isl). The interaction 
of C with the (/c + l)th atomic sample begins as soon as that with the kth sample ends. Successive 
atoms are thus separated by the time interval T. We denote by pk the cavity state just after 
interacting with the kth atom and tracing over its irrelevant final state. The field density operator 
Pk is thus given by: 

Pk = M^^pk.iM^^ t + My-pfc.iM^- t . (6) 

We aim to stabilize a pure pointer state poo = |'0oo)('0oo|7 which must be a fixed point of this 
Kraus map. The right-hand side of Eq. (p|, with pk-i = Poo, is a statistical mixture of two pure 
states. It is then a pure state only if its two terms are proportional to each other. Thus, |'0oo) must 
be an eigenstate of both M^^ and M^^ . State stabilization by reservoir engineering amounts to 
tailoring a Kraus map for the field from a constrained physical setting. 

We have shown in [32 that it is possible to engineer the atom-cavity interaction so that the 
Kraus map leaves invariant the states ~ g-^^^H^ |q,^^ [^ which a and ^k ^k in Eq.Q can be chosen 
at will. Explicitly, we build Ut by sandwiching a resonant interaction ((5 = for t G [— 1^/2, tr/2]) 
symmetrically between two dispersive interactions with opposite detuning: 5 = 5q before the 
resonant interaction, 5 = —Sq thereafter. We choose a positive So value for the sake of definiteness. 
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Figure 3: Time profile of 5 (difference between the tunable frequency cjq of the atomic transition 
and the fixed frequency ujc of the single cavity mode; bottom frame, full line) and Vt{vt) (bottom 
frame, dashed line) during cavity crossing by one atomic sample; we take t = when the atom is 
at cavity center. Top frame: scheme of the propagators corresponding to the successive steps in 
the composite interaction. 



This timing is illustrated on Fig. [3J Each resonant or dispersive interaction phase is charac- 
terized by a set of parameters that we denote q = (ti, ^2, '^, (^o) where ti is start time, ^2 is stop 
time. The corresponding propagators are denoted Ug. 

To compute these evolution operators, we use the fact that each two-dimensional subspace 
spanned by (|^, n -F 1), |e, n)) is invariant under the action of Hjc- The state \g, 0) does not take 
part in the evolution. We can thus view the action of \Jq as photon-number-dependent rotations 
acting on a set of Bloch spheres Bn {n = 0, 1, ...), with |^, n + 1) at the south-pole and |e, n) at the 
north-pole defining their Z-axes. These rotations can be decomposed as general rotations around 
the X-, Y- and Z-axes of the Bloch spheres. We use the notation /n = /(N) = ^^ f{n) \n){n\ = 
Zln fn l^)(^l ^o^ arbitrary functions / of n, with the fundamental property 

a/(N)=/(N + I)a. (7) 

We can then cast these rotations as: 

X(/n) = |f/)(5|cos(/N/2) + |e)(e|cos(/N+i/2) (8) 



Y(/n) = |5){5|cos(/N/2) + |e)(e|cos(/N+i/2) (9) 



Z(/n) = \g){g\e'f^/^ + \e){e\e-'f-^^/\ (10) 

As shown in Sections b^ and Ul Y(/n) with f{n) proportional to y/n -\- 1 corresponds to a resonant 
interaction and Z(/n) with /(n) proportional to n corresponds to a non-resonant interaction in 
the dispersive regime S :^ ft. See Appendix [A| for more details. 

3 Engineered reservoir for coherent state stabilization 

The coherent state \a) is obtained by the action of the displacement operator Dc^ _ g^a -a a 
onto the vacuum i4|: 

|a)=D«|0) =e-H/2^^|n). (11) 



We show here how a short resonant interaction (^ = 0) with weakly excited atoms provides an 
engineered reservoir, whose pointer state is close to a coherent state [32 . 

Stabilization of coherent states is not an amazing feat. A classical radiation source weakly 
coupled to the cavity directly generates a coherent state. This is a routine operation in microwave 
CQED experiments. However, the situation described in this Section is an essential building block 
for the stabilization of more complex nonclassical states, as shown below. Moreover, it is an 
interesting micromaser situation ^81 136] , in which the small excitation of the atoms leads to a 
finite energy in the steady state even though the cavity is assumed to be lossless. 

We consider a resonant interaction over a time interval t^, corresponding to the parameter set 
r = (— 1^/2, t^/2, 'u, 0). Following Appendix \A\ the associated propagator is: 

U, = Y{e'^) (12) 

= |ff)(5|cos(^N/2) + |e)(e|cos(^j;j+i/2) 

with 



rU/2 

= OrVn, Or= n{vt) dt . (13) 

J-tr./2 



This readily yields: 



M^" = cos(f ) cos(ggj/2) + sin(f ) '"'"^^/^^ a^ 



MY- = sin(f) cos(^?,+i/2) -cos(f)a ^'"^J|/^^ . (14) 



A pointer state of this resonant reservoir must be an eigenstate of both M^ and Mg given by 



Eq.(14). Let us expand it over the Fock states basis, IV^oo) = X^ V^nl^)- We get a condition on the 



coefficients ipru for n = 0, 1, 2, ...: 

sin(^;+i/2) V^^+i = tan f (1 + cos(^;+i/2) ) ^^ • (15) 

This relation allows to compute all the ipn starting from ipo ^ 0, except if sin(^J^/2) = for some 
m. This condition corresponds to the existence of a trapping state |m — 1) [37 , which is then 
uncoupled from |m). The Hilbert subspaces corresponding to the photon numbers < (m — 1) and 
to those > m are then decoupled, such that the steady state depends on the initial conditions. 
Since all states considered in the remainder of the paper have an energy lower than 20 photons, 
we arbitrarily truncate the Hilbert space to n < nmax = 50. We can thus avoid trapping states 



by choosing small dr values such that sin(6>JJ^^;L/2) 7^ for all < n < nmax- Dividing (15) by 
sin(^J^^;L/2) then gives the recurrence: 

, tan(2i/2) 

^"+^ = tan(^;^,/4) ^" ' (^*^) 

which defines a unique normalized pointer state. 

For (^nmax/^)^ ^ -'-' ^^ recurrence (^16| reduces to V^n+i ^ an(n/ J ^^^ which defines a 

coherent state \oloo) with Ofoo = — ^ ^ (compare with the last member of Eq. (pT|). 

This value of a^o can be retrieved by a simplified reasoning as in [32] . Assume that the cavity 
already contains a large coherent field of amplitude a ^ 1. The incoming atoms then undergo 
a resonant Rabi rotation in this field, with an atomic Bloch vector starting initially towards the 
south pole of the Bloch sphere. The Bloch vector rotates by an angle — 6>^a, such that if drOt < 2u 
(resp. OrOL > 2u) the final atomic state has a lower (resp. larger) energy than the initial one, i.e. in 
average gives energy to (resp. draws energy from) the field. This energy exchange thus stabilizes 
a field with amplitude Ofoo = ^u/Or- 

We have numerically examined the fidelity F = |('0oo|<^*)P of the pointer state |V^oo) defined 



by Eq. (16) with respect to a coherent state \a^) of the same mean photon number \a^\ = 
(V^oolNl'^oo)- Figure Qa) represents that mean photon number in gray scale, for a range of 
parameters u^Or delimited such that the fidelity F is larger than 99%. We limit the plot to 
6r < (27r)/y'nmax ^ 0.88 to avoid trapping states, and to Or > 5 tan(ii/2)/y'nmax to remain in 
the truncated Hilbert space (top left corner cut off). The coherent state approximation for iV^oo) 
remarkably holds for a range of u and 0^ much larger than that predictable from the qualitative 
discussion above. 

Convergence towards |aoo) can be simply analyzed in the limit of small u^Or. Expansion of 



Eq. (14) to second order in u^ 6^ yields the Kraus map: 



uO 

Pk+i ^ pfc + ^ ( [a^ Pk] - [a, pk] ) (17) 

02 

-^ (Npfc + P/cN - 2ap/ea'^) . 
o 

It can be simplified by letting p = D_q,^ P^c^oa such that p = |aoo)(<^oo| corresponds to p = 



A few calculations show that Eq. (17) transforms into 



02 

Pk+i = Pk- Y ^^Pk + P/cN - 2a.pka.^) . (18) 



This is a finite difference version of the standard Lindlblad equation, describing the damping of 
an harmonic oscillator coupled to a zero temperature bath. It efficiently drives any initial state 
towards the vacuum, with a relaxation rate proportional to 0'^. This analogy shows that the initial 
map [Eq. ( p!7| )] drives any initial cavity state towards a coherent state |aoo) with aoo = 2u/6r. 
Smaller Or values, i.e. shorter interaction times of each atom with the field lead to a higher energy 
pointer state (for a given u)^ but to a lower convergence rate (independently of u). 



Similar conclusions are directly reached from Eq. (17) by assuming that the field is at any 
stage during its evolution towards equilibrium in a coherent state with amplitude a^. Using 
simple second-order approximations in u^Or [38 , we find that this amplitude evolves as: 

c^fc+i = (l-^'/8)afc + ^x^./4. (19) 

This first-order system has the explicit solution a/e = (1 — 6>^/8)^ (ao — <^oo) +Qfoo starting from a^ 
at A: = 0. Noting that log |(aoo|<^/c)P = — |<^/c — <^ooP, the fidelity indicator log | log |(aoo|<^/c)Pl = 
log |ao — Q^ooP — Aconv k decreases linearly in k towards — oo. The slope Aconv = 2| log(l — ^^/8)| 
measures the exponential convergence speed of \ak — Q^ooP, which increases with Or and is inde- 
pendent of u. 

Numerical simulations of Eq. (|6| with the exact Kraus map [Eq. ([l4|] vindicate this approxi- 
mate analysis. Figure[4]^b) shows the evolution of log | log('0oo|p/c|'0oo)| as a function of the number 
of atom-field interactions /c, starting from the vacuum po = |0)(0|7 with the real Kraus map as- 
sociated to Ur- The evolution is linear, as predicted by the simplified model. We have checked 
that this linearity holds for a large range of parameter values: it is only for large Or^ that the 
curve bends slightly upwards for the first few /cs. This allows us to use the slope Aconv of that 
approximate line as a measure of convergence speed. Fig. [4]^c) shows the dependency of Aconv in 
Or^ for two different u values: i^ = 0.1 (dotted curve) and u = 1 [dashed curve, which does not 
extend to low Or values, according to the accessible domain on Fig. [4Fa)]. They closely follow 
the simplified model (full line) , which is independent of u and slightly overestimates convergence 
speed. 

4 Kerr Hamiltonian simulation in the dispersive regime 

We now discuss the case of a full composite interaction, with the detuning profile 5 it) represented 
on Fig. |3] The full propagator 

Ut = U,, U, U,, , (20) 

is the concatenation of three unitary operators corresponding first to the dispersive interaction with 
di = (— T/2, —tr/2, V, So) as the atom enters the cavity, then to the resonant one with q = r around 
cavity center, then to the last dispersive interaction with d2 = {tr/'^^ ^/2, v^—5q) as the atom leaves 
the cavity. The exact expression of U^ is given by Eq. ([l2| . The dispersive propagators XJ^i and 
\Jd2 are computed in Appendix [Aj assuming that 5q and v satisfy the adiabatic approximation 
condition [Eq. (|46])]. This condition merely expresses that the interaction Hamiltonian Hjc varies 
slowly (through Vt{vt)) in comparison to the differences between its eigenfrequencies. Thus, each 
initial eigenstate of Hjc remains an eigenstate and only accumulates a phase proportional to its 
eigenenergy. 

To get a simple insight in the physics of the stabilization of nonclassical states, the present 
Section focuses (like [32j) on the simple case in which the two nonresonant interactions take place 
in the dispersive regime, i.e. 5q^ VLq. This avoids spurious population transfers in the dispersive 
phase, as atomic levels dressed by the cavity field almost coincide with the bare levels |e) and \g). 
The dispersive propagators, deduced from Eq. (52) in Appendix [A} then write: 



it ^ ^{J.d 



Urf, « m « z(</.^) (21) 



with 0N = 0^N ■ 
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Figure 4: Reservoir with interaction U^. (a) Mean photon number ('0oo|N|'0oo) of the pointer 
state IV^oo)- Grayscale axis is hnear in \/ {il^oo\^\'4^oo)' The shaded zone is dehmited such that 
the corresponding states have at least a 99% fidelity KV^ool^^*)!^ to a coherent state \a^) of same 
mean photon number |a*p = ('0oo|N|?/;oo). On the top left corner, pointer states have significant 
population outside the truncated Hilbert space. On the top right part, KV^ooIci^*)!^ drops to ~ 98% 
as u approaches 7r/2. (b) Evolution of the fidelity indicator log | \og{iljoo\pk\'^oo)\ as a function of the 
number of atom-field interactions (i.e. Kraus map iterations) /c, starting from vacuum po = |0)(0|- 
We have arbitrarily set u = 0.5 and 6r = 0.4, for which ('0oo|N|?/;oo) = 6.21. (c) Convergence rate 
Aconv as a function of 6^^ for ix = 0.1 (dotted curve) and u = 1 (dashed curve). Dependency in u 
is small. We also represent (full curve) the analytic result of the simplified model [Eq. (Il9|]. This 
model slightly overestimates the convergence speed. 
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where 0-y = 1/(2 (5o) S_x L ^'^{i^t)dt is a phase shift per photon and (f>(^ = So{T — tr)/2 reflects 
the free atom evolution in the interaction representation at cavity frequency. 
The full propagator then writes 

Ur«Ud = Z(-0 U, Z(^^) 

= \g){g\ cos e^^/2+\e){e\ cos e-^+j2 
sin 6*^/2 ^i(^^]sr+0^) 



|e){ff|a \£^ e 



-|5){e|^^^e-W^N+,c)at, (22) 



where 6^ is defined by Eq. (13). The opposite dispersive interactions have no net effect when 
the atom remains in the same state during the resonant interaction 11^. In contrast, Z does not 
commute with terms in which the atomic level changes in the resonant phase. For these terms, 
the dispersive phase shifts add up. The global evolution XJ^ thus associates a phase shift to each 
term of Ur that changes the field energy. An increasing field energy corresponds to a decrease 
of the field phase and vice versa. These correlated phase and amplitude shifts suggest that \Jd 
might stabilize a coherent state distorted by amplitude-dependent phase shifts, a situation similar 
to that encountered during the propagation through a Kerr medium. 

The action of the atom can indeed be expressed by an operator acting on the field only. Let 
us define the Hermitian operator h^ by: 

h^ = (/)^(N' + N)/2 + ^^N . (23) 

With the commutation identity [Eq. ([7|] we have e~*^N ae*^N = a 6*'^'^^^+'^^^ and 

Vd = e"'^^ U^ e'^^ . (24) 

Thus, JJd is equivalent to Ur modulo a basis change on field state alone defined by the unitary 
operator e~'^^^ . 

In other words, when p evolves under the Kraus map associated to (M^^ , M^^^), p^ = 
e*^N p g-*^N ^ evolves under the Kraus map associated to (M^'^ , My^^). It follows from Section 3 
that p^ converges towards a coherent pointer state |Qfoo). Therefore, p converges with the same 
convergence rate towards a nonclassical pointer state ex.-p[—ih^] \aoo). 

The effective Hamiltonian h^/tx is equal to the Kerr Hamiltonian H^, with ^k^k = ^7/2 
and Ck^k = (</>c + 4^i/'^)' The engineered reservoir thus stabilizes the nonclassical pointer states 
^-itK^K \aoo) which would be produced by propagation through a Kerr medium (see Fig. 11]). 
Tuning T and 5q allows us to choose (j)^ at will. We can thus prepare and stabilize a whole 
class of such states, as described in Section fl] In particular, for (j)^ = tt, we get the MFSS 
kttoo) = il^oo) + ^ |-<5^oo))/a/2 with aoo = e~^^^'^^^'^^ aoo- Note that the stabilization of this two- 
component MFSS is the most demanding experimentally, since it requires the longest dispersive 
interaction time. 

The discussions in this Section only apply in the limit of small VL/Sq. Reaching notable (/)^ in 
this case requires a large dispersive interaction time (T — t^)/2, that can be prohibitive. First, 
in the experimental context of Fig. [2I T = ?tw /v is limited by the minimal achievable atomic 
velocity (a few tens of m/s in the ENS setup). Second, a larger T means less frequent atom-field 
interaction and thus a weaker reservoir, implying a less efficient protection of the target state 
against decoherence induced by cavity relaxation. 

5 Regime of arbitrary detunings 

We thus discuss now dispersive interaction with moderate VL/Sq values, which allows to reach 
significant dispersive effects within moderate interaction times. We therefore use a more precise 
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expression of the propagator for the nonresonant interactions (parameter sets di and ^2), by 
applying the adiabatic approximation to the actual dressed states (instead of \g^n + 1) and |e,n) 
as in Eq. (21) when assuming Sq ^ Qq). Developments detailed in Appendix [A| lead to: 



Ut « Ue = Z(-(/.n) X(^n) Y{0^^) X(^n) Z(,^n) , (25) 

with 

r-U/2 

^n = So ^l^n(Sl{vt)/5oYdt , (26) 

J-T/2 

tane„ = iMV« withe„G(^,f). (27) 

We recognize in this expression the central resonant interaction evolution operator, Y(6>J^), and 
the two phase-shift operations accumulated during the non-resonant interactions (Z(— (/)n) and 
Z((/)n))- Note that here, unlike in Section [4J (j)ri is a nonlinear function of n. The remaining 
two X(^n) operators reflect the fact that the atomic energy eigenstates do not coincide with the 
dressed levels at ±t^/2, when the atomic frequency is suddenly switched. Note that we neglect 
two similar transformations which appear in principle when the atom gets first coupled to the 
mode and finally decoupled from it, since the atom- field coupling is then quite negligible. 
Some tedious but simple computations exploiting Eq. ^ allow us to write: 

Ue = |<7)(5|cos(^^/2) + |e)(e|cos(^^+i/2) 



Here, 6% G [0, 27r) is defined by 



■ |g)(e|e-»<^N sin(y) ^t, (28) 



N 



cos(0y2)=cos(0y2) cosCn. (29) 

Introducing ^S] 

X'n = angle[ sin(^y2) - i cos(^y 2) sin e„ ] , (30) 

we define the composite phase as (j)^ = (j)-^ + Xn- 



Comparing Eqs. (28) and (22), we finally get: 

U, = Z(-0^) Y(0^) Z{4>h) . (31) 

This expression of Uc has the same general form as that used in the dispersive case (Section [4|. 
Angles ^^, (f)^ replace ^J^, (j)^ respectively. We now show that with these adaptations, most of the 
conclusions of the previous Sections still hold. The reservoir in realistic situations indeed stabilizes 
the nonclassical states \iIj) ^ e~*^^^^ \a). 

5.1 Effects of Y(e^) and Z(±(/)^) 

Let us first consider a reservoir of atoms whose interaction with the cavity would be described by 
Y{0^). Note that this situation is not physical: the Y{6^) evolution operator is no more than 
a convenient mathematical factor appearing in the expression of the complete evolution operator 

In analogy with Section pj the pointer state \^|Joo) = ^ V^nl^) corresponding to this fictitious 
interaction is defined by the recurrence relation: 



12 



for n = 0,1,2,... . Equation ([29| ensures | cos((9^+i/2)| < 1 Vn, therefore < 6>^+i/4 < 7r/2. 



Moreover lim^^+oo ^n ~ '^^ such that Eq. (32) always yields a well-defined finite energy state as 



soon as |tan(ii/2)| < 1 i.e. \u\ < 7r/2. For large n values, the recurrence (32) is approximated 
by V^n+i = tSiU^ ipri and ipn converges exponentially towards with ^^ n?/^^finite. The energy 
exchange resulting from the fast commutation of the atomic frequency near the cavity center thus 
removes the possibility of trapping states. 



Note that even in the absence of the central resonant interaction, with ^J^ = in Eq. (29) 



relation (32) defines a unique pointer state with finite energy. It is thus in principle possible to 



simplify our scheme by using only two dispersive phases with opposite detunings. 



For Or and Qr/So smah, 2tan(l9^/4) ^ 0^/2 ^ 6>cV^/2 with 0^ = Joi^(^Y. The 
pointer state is thus close to a coherent state |aoo), as in Section [sl with the amplitude Ofoo = 
4 tan(ii/2) /6>c. Convergence arguments similar to those of Section^ (effective Lindblad master 
equation) can be given. The convergence rate is now proportional to O^. We conjecture that this 
convergence is valid for any u with {) <u < 7r/2, ^^ > 0, l^r > and 5q > 0. 

Figure [5] presents numerical estimations of the field pointer state \^|J 00) stabilized by a hypotet- 
ical reservoir using interaction Y{0^). For all represented parameter values, fidelity KV^oolo^*)!^ 
to a coherent state \a^) of the same mean photon number (|a*p = ('0oo|N|'0oo)), is at least 99%. 
That mean photon number is represented by the grayscale. The larger value ftr/^o = 1/2.2 used 
on Fig. [sja), does not allow to reach as high mean photon numbers as the value Qr/^o = 1/10 of 
Fig. [5|b). Small ftr/^o however are more subject to undesired population of the high- lying Fock 
states, reminiscent of the trapping states, for large Or and u. This explains the smaller domain 
where fidelity is larger than 99%. With the particular values used for Fig. [ije) (black dot on 
Fig. [sja), with Qr/^o = 1/2.2, = 7r/2 and u = 0.457r), fidelity to a coherent state is almost 
99.9% and (V^oo|N|V^oo) = 2.96. 

We now examine the infiuence of the Z(±^^) operators on the pointer state defined by Y(^^). 
A first observation is that it does not modify the photon number populations, since it commutes 
with N. Thus, the energy of the field pointer state, for a reservoir with composite interaction, 
is entirely determined by the parameters in Y(6>^), as represented for example on Fig. Is] Let us 
define the Hermitian operator h^ by the recurrence relation: 

/in+i -hl = C+1 , (33) 
for n = 0, 1, 2, ..., with an arbitrary Hq. Using Eq. (JTl as in Section [4] yields 

Uc = e-^^N Y{0^) e^^N . (34) 

The pointer states of Uc are those of Y(^^) transformed by the unitary e~*^N. The choice of 



Hq for solving Eq. (33) is thus physically irrelevant, as it corresponds in Eq. (34) to two constant 
opposite phases that cancel out. The operator h^ here plays exactly the role of h^ in Sectional 



The only difference is that, as (j)^ is nonlinear, h^ is defined with the discrete integral (33). If < 
is nearly linear in n over the relevant photon numbers [dominant photon numbers in the pointer 
state iV^oo) associated to Y(6>^)], then h^ is nearly quadratic and the situation of Section [5] is 
recovered. The reservoir stabilizes nonclassical pointer states IV^^) = e~*^N 1?/;^) ^ e~*^^^^|a) 
with tK depending on the parameters governing (j)^. 

5.2 Choice of the reservoir operating point 

We now use this detailed description of the reservoir to justify the choice of operating parameters 
leading to the generation of the two component MESS presented in figure file): u = 0.45 tt. 
Or = 7r/2, 'u = 70 m/s, S = 2.2(^0- Note that the state in Fig. file), with ^ 2.7 photons on the 
average has been computed with a finite cavity lifetime Tc = 65 ms. The same computation in an 
ideal cavity leads to an average photon number equal to 2.96, see Fig.|5| The two- component MESS 
corresponds to the largest effect of the dispersive interaction, and hence to the most demanding 
experimental conditions. 
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Figure 5: Mean photon number {^oo|N|^oo) of the pointer state \tpoo) stabihzed by Y{0^), •with 
5o = 2.2 fir (a) and i5o = 10 1^^ (b). Grayscale axis is hnear in \/{ipoo\^\ipoo)- The 2-component 
MFSS on Fig. IT] uses u = 0.97r/2 and 9r = tt/2 •with Sq = 2.2 fir, corresponding to the black dot 
on (a), for •which {^oo|N|-i/'oo) = 2.96. The shaded zone is delimited such that all corresponding 
states have at least 99% fidelity |(^oo|q;*)|^ to a coherent state |a«) of same mean photon number 
(|a*|' = (^oo|N|^oo)). 
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Figure 6: (a) Kerr-effect-inducing Dcj)^ = 4>^^i — 0^ as a function of photon number n. Since (j)^ 
is undefined, we start with D(l)\. We set Or = 7t/2. (b) Corresponding velocities v: for each (^o/^r, 
we adjust v to have Dcj)^ = tt at n = 2.96. That value is chosen to cover the parameter values 
V = 70 m/s, So/^r = 1/2.2 (black dots) used for the 2-component MFSS in Fig. [l] An ideal /i^, 
proportional to H^, requires D(j)^ constant in n. Small So/Qr values are disadvantageous for this 
criterion, but allow higher velocity and hence more frequent reservoir atoms for a same mean Dcj)'^. 



The chosen parameters are the result of a tradeoff between contradictory requirements. First, 
the composite phase shift 0^ must be nearly linear in n over the useful photon number range, 
with a slope of n per photon. Second, the time of convergence towards the steady state needs 
to be much shorter than the decoherence time (Tc/5.6) of the target state due to unavoidable 
cavity relaxation. Linearity of (/)^ improves with larger So/Qr- The tt phase shift per photon 
condition then requires very long atom-cavity interaction time, in clear contradiction with the 
second requirement. 

The tradeoff is further examined on Fig.|6] Figure |6ja) evaluates the linearity of (j)^ by showing 



D^n 



^n+1 



^n+1 



2h^ for different parameter values. Once again Or has little 



influence and we set it to 7r/2. For each value of ftr/6o^ we adjust v to have D(j)^ = tt for the same 
mean photon number n = 2.96 (by interpolation). That value is chosen to cover the parameter 
values of Fig. me), represented by a black dot on Figs.plandpJ As expected, Dcj)^ is quite constant 
for moderate photon numbers in the dispersive region 60/ Qr ^ 1- This corresponds however to 
unrealistically small atomic velocities, represented on Fig.[6]^b). In the region of low So/Qr values, 
a D(j)^ ^ TT at n = 2.96 can be reached with larger v, but Dcj)^ varies more rapidly with n. This 
variation is nevertheless sufficiently weak in the range 2 < n < 5 for So/Qr ^2.2, corresponding 
to the 'U = 70 m/s that is used for Fig. [lie). 

Let us now examine the overall reservoir fidelity and the convergence rate Xconv from the 
vacuum towards the target state, as defined in Fig. |4l We choose as free reservoir parameters Or 
and So/Qr- This choice sets the value of v [see Fig. Jwb)]. Then u is adjusted so that the target 
mean photon number is 2.96 (see Fig. Isja)). Figure IWa) shows the ratio Xconv /T, where T is the 
total interaction time of each atom with the cavity. This ratio is the real convergence rate in s~^ 
units. For Fig. [ije), since the expected target state decoherence time is of order 65/5.6 ^ 10 ms, 
we choose a parameter set 60 /fir = 2.2 and Or = 7r/2, corresponding to a 1400 s~^ convergence 
rate (black dot on Fig. JTJ^a)). The fidelity with respect to an ideal MFSS with the same energy 
is shown on FigJTFb). Our choice of parameters does not correspond to a maximum fidelity due 
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Figure 7: (a): Convergence rate Aconv/^ giving the slope, in time units (s~^), of the convergence to- 
wards the reservoir pointer state IV^^), according to log | \og{ip'^\pk\il^^)\ = log | \og{ip^\po\tlj^)\ — 
Aconv ^ (see Fig. [4|. For each Or and l^r/^o, we adjust v as in Fig. |6] to keep Dcj)^ ^ tt, 
and 1/ to keep (V^^INI?/^^) = 2.96; this reference comes from the values Or = 7r/2, i/ = 0.457r, 
Qr/^o = 2.2, '^ = 70 m/s (black dot) used for Fig. [ije). Time T = 3w/v between consecu- 
tive atoms changes as we adjust v. (b): Fidelity of the same \^|^'^) to a 2-component MFSS 
Wa ) — il^oo) + ^e*^ \-aoo))/V^j where we tune aoo and < /3 < tt to optimize fidelity. It turns 
out that |/3| < 0.005 for most parameter values, while lofooP decreases as fidelity decreases, below 
2.7 for the lowest values of Qr/^o- The black dot marks the case of the 2-component MFSS in 
Fig. [l] For Or values larger than those represented, no u value stabilizes a mean photon number 
2.96, see also Fig. [5j The two plots together illustrate a tradeoff between fidelity in absence of 
decoherence and convergence speed. 

to the variation of Dcj)^ in the useful n value range (see Fig. [6|a)). However, we get an excellent 
95% fidelity. 

Figures [S] left and right respectively present the Wigner functions of the steady state MFSS 
obtained with this parameter choice, and of a theoretical superposition of two coherent states with 
opposite phases and same total energy. The slight distortions of the quasi coherent components 
in the pointer state MFSS are due to the non-linearity of the phase shift (j)^. 



6 Decoherence and experimental imperfections 

The choice of operating parameters performed above has been based on a rough estimate of the 
action of decoherence. We now show how the reservoir allows us to stabilize MFSS with a high 



fidelity, in presence of cavity relaxation due to a zero-temperature environment (Section 6.1). In 
Section |6.2[ we study the robustness of the scheme against other experimental imperfections by 
numerical simulations. 
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Figure 8: Wigner function illustrating a stabilized 2-component MFSS (colorbar as in Fig. fTl). 
Left: |'0^)('^^|, pointer state of our reservoir with composite interaction. Parameter values are 
6r = 7r/2, u = 0.457r, 8q = 2.21^0, '^ = 70 m/s, i.e. those used for figure [ije), except Tc set to 
infinity here. Right: target state |cq,)(cq,|. 



6.1 Cavity relaxation 



First, we consider the field evolution with a simplified model. It is obtained from equation (19) 
for a coherent evolution, sandwiched between two dispersive transformations (Eq. ([21])): 



Pk 



(1 - el/d>)ak 

\(y-k){(^k\ 



t/2 N^ 



f^^^ 



uOr/^ 



r/2 N^ 



(35) 



In a Monte-Carlo approach, the evolution of the field density matrix due to cavity relaxation 
can be represented as a succession of quantum jumps described by the annihilation operator a, 
occurring at random times and interrupting a non-unitary deterministic evolution of the field state 

m- 

The action of a on an MFSS |ca) writes: 

a|cc,) (X \c_oc) ' 

Since the loss of photons in the environment cannot be detected, an initial |cq,) state rapidly evolves 
in the absence of the reservoir into a statistical mixture of \ccy) and |c_q,), i.e. into a mixture of 
\a) and |-a). When the reservoir is present, it drives \c.cx) back to \coc) after each jump. If the 
reservoir-induced convergence time is much shorter than the average interval between two jumps, 
then the field is mostly close to |cc^). 

This simple description suggests to seek a solution for the steady state with decoherence under 
the form p'^ = J p{z)\z){z\ dz. This is a statistical mixture of coherent states \z) with real 
amplitudes z weighted by /i(^). 

In the absence of cavity relaxation, the evolution of p'^ in the simplified model [Eq. ( |35| )] can 
be viewed as a discretization of the Lindblad master equation: 



^-'^ = [/3at 



dt 



/3ta, p'^] - f (Np^^ + p'^N - 2ap'^at) , 



(36) 



with (3 dt = uOr/4: and K,dt = 6>^/4. Eq. ([36| describes the evolution of the field mode coupled 
with a classical source with an amplitude /3 and damped at a rate k. At long times, p' converges 
towards a coherent state |aoo) with aoo = 2/3/ai:, see e.g. [4^. Note that in Eq.(36), the damping 
rate k, is induced by the atomic reservoir and not by cavity relaxation. Since //^follows (36), p' 
follows 



dti^ 



[ate- 



tN 



^^a, p'] 



f (Np' + p'N - 2e'"^ap'ate-^"^) , 
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where we can assume, up to a change of phase, that f3 is real and positive. 

We now add to this simple model a thermal environment at zero temperature that induces 
decoherence of the field with the cavity lifetime Tc = l/t^c- This adds the usual Lindblad terms 
to the second member of the previous equation and p' now obeys: 



f^p' = p[a^e-'^^ - e^^^a, p'] 



- f (Np' + p'N - 2e*"^ap'ate-'"^) 

- fCNp' + p'N - 2ap'at) . (37) 

In the Kerr representation, p^^ then evolves according to: 



^-'^ = /3[at-a, p'^] 



dt 



^{Np'^ + p'^N - 2ap'^a^) 

- /^c(ap'^at - e'^^ap'^R^e-'^^) . (38) 



Without the terms in the third line of Eq. ( |38| ), we would get Eq. (|36| with hz replaced by hz-\- hZc- 
This would yield a coherent steady state of amplitude a^ = Qfoo/(l + v) with r] = (4T)/(^^Tc). 



The whole equation (38) leaves invariant the set of mixtures of coherent states with real amplitudes 



in [— Q^^, <^^]. We therefore search for the stationary solution under the form: 

P'^= r°° H{z)\z){z\dz. (39) 

As shown in the Appendix [B} this yields a solution: 

Kz) = A^o ^ ^^-3^ ^— , (40) 

with Tc = 2hZc/{i^ + ^c)- The normalization factor /io > ensures that J^^^ p{z)dz = 1. In any 
case, p{—a^) = 0. For small hZc^ we have lim^^Q^c^ p{z) = +oo and p^^ is close to the coherent 
state \a^). For large hZc^ cif^ tends to zero and thus the field steady-state becomes close to the 
vacuum. 

We now compare this simplified model to the actual reservoir in the presence of relaxation. 
Figure [9] illustrates the reservoir-induced convergence after a quantum jump. The leftmost column 
shows the Wigner function of p' during this recovery process for the simplified model ( |35| . We 
start as \ca) (upper left frame). Immediately after a jump (second frame in the leftmost column), 
the state is |c_q,). Successive snapshots of the recovery procedure are presented in the next frames. 
We neglect here the action of cavity relaxation during this recovery process. Note that after ^ 4 
reservoir atoms, the state is the vacuum, from which \ca) is gradually recovered. 

The second column depicts the evolution of p'^. In this representation, the initial state is the 
coherent state \a) (first frame). It jumps to \-a) (second frame), and then gradually evolves back 
towards \a) according to Eq. ([l9|, staying coherent at all time. 

On the third column, we show the Wigner functions of the actual cavity state p induced by 
our reservoir, whose dynamics is governed by the Kraus map associated to Uc- The last column 
shows the evolution of p^ = e*^N pe~*^N. We observe that p^ and p'^ follow qualitatively the 
same path. The main difference is a notable distortion of p^ when the field amplitude is near zero. 



In figure [lOj we plot the two marginal distributions of the Wigner functions for p^^ and p^ 
along the real and complex quadratures. The reservoir steady states p'^ and p^ approximately 
correspond to the quantum Monte Carlo average of the trajectories depicted in Fig. [9J Figure [lO] 
features dominant peaks which suggest that the field is mostly close to the target. The distortions 
with respect to a coherent state visible on the fourth column of Fig. |9J lead to a plateau or bump 



on the marginal distributions of p^. We nevertheless observe that our simplified model [Eq.(37)] 
captures the main features of the infiuence of decoherence. 
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Figure 9: Evolution of the cavity field coupled to a reservoir stabilizing a 2-component MFSS, 
immediately after a relaxation-induced photon loss. Parameter values are those used for Fig. [8J 
The frames are labelled by the number of atomic interactions. A photon loss out of the reservoir 
pointer state occurs between frames labelled 0~ and 0+. Left two columns: simplified model, 



described by Eq. (35). We show the Wigner functions of both the cavity state, p' ^ on column 1 
and of p'^ on column 2. Right two columns: same plots for the actual reservoir characterized by 
Uc (p on column 3 and p^ on column 4). 
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Figure 10: Steady state of the cavity field coupled to the atomic reservoir and to a relaxation- 
inducing environment with Tc =65 ms. Top [resp: bottom]: marginal distribution of the Wigner 
function along the imaginary [resp: real] quadrature for the simplified model (dashed line) and 
for our reservoir (solid line). These states correspond to the quantum Monte Carlo average of the 
sequence presented in Fig. [9J 
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6.2 Experimental uncertainties 

We have performed extensive numerical simulations to assess the robustness of the reservoir versus 
uncertainties in the experimental settings. The evolution operators during the interaction are 
computed exactly from the Hamiltonian Hjc [Eq. (|3|], using the quantum optics package for 
MATLAB [41]. Hilbert space is truncated to the 60 first Fock states. We take into account 
separately the atom-cavity coupling ruled by Hjc and the relaxation of the cavity mode (modeled 
in the standard Lindblad form). This simplifying approach holds since T <^Tc. 

We use as a reference the generation of a two-component MFSS containing 2.7 photons on 
the average (thermal environment with Tc = 65 ms and a mean number Ut = 0.05 of blackbody 
photons per mode, S = 2.21^0, t' = 70 m/s, u = 0,457r, tr = ^ /is, see Fig. [ije)). We take into 
account the randomness of the Rydberg state preparation [32]. In each atomic sample, we excite a 
random number of atoms obeying a Poisson law with an average Pat- We take a low value Pat = 0.3 
such that, in a first approximation, we only get samples containing zero or one atom (we examine 
later in this Section the influence of samples containing two atoms). Note that these are the 
conditions used for figure [ije). In each case, we compute, for a slight change in the experimental 
settings, the variation of the fidelity of the pointer state w.r.t. an ideal optimized two-component 
MFSS. For the reference set of parameters, the steady-state fidelity is 70%. 

For velocities in the 66 < v < 74 m/s interval, the fidelity is only slightly altered, varying from 
65% to 70%. It is thus insensitive to a velocity dispersion in the 10% range, well below the values 
achieved in the experiment. 

The fidelity is also quite insensitive to a slight mismatch in the values of the detuning for the 
two dispersive interactions. Assuming that 5 takes the value ai x 2.2(^0 in the first dispersive 
period and —a2 x 2.2 Qq in the second, the fidelity drops by at most 10% when ai and a2 vary by 
up to 10%. The latter cover far more than the actual uncertainty on the atomic frequency. 

We have slightly offset the timing of the resonant interaction, shifting it in time by 5t and 
keeping tr constant. A shift of up to 1 /is (well above the 10 ns timing accuracy) has no effect 
on the fidelity at the 1% level. The fidelity is also quite insensitive to a finite rise time for the 
voltage controlling the atomic Stark effect in the cavity, and hence to a finite commutation time 
for the detuning S. Using an exponential relaxation model, and setting tr to maintain a constant 
Or value, we find that the fidelity is unchanged for commutation times up to 200 ns, in the range 
of accessible values. 

We have also studied the sensitivity to the atomic samples containing two atoms at the same 
time. We decide randomly for each sample the actual number of atoms, A/'^, according to a 
Poisson distribution with the average value Pat^ truncated above Na = 2. For two- atom samples, 
we integrate the exact equations of motion, assuming an identical coupling of both atoms to the 
mode. This condition is realized in the experiment, since the maximum separation between the 
atoms in a sample is, below 1 mm in C, much smaller than the wavelength - 6 mm - or than the 
mode waist w. 

We observe that the two-atom events do have an impact on the fidelity. For Pat = 0.3, the 
energy of the prepared cat decreases down to 2.4 photons on the average and the fidelity is reduced 
to 66%. For larger Pat values, the decrease is more important and the fidelity reduces to 34% for 
Pat =0.5 (for larger pat^ the simulation should also include 3- atom samples). 

If we consider an unrealistic reservoir involving samples containing always two atoms, we get 
as steady-state a large two-component cat, with 4.8 photons on the average and a fidelity of 65%. 
In the real situation, this two-atom engineered reservoir interferes destructively with the operation 
of the one- atom samples, leading to reduced average energy and fidelity. 

When we reduce Pat below 0.3, the fidelity and the energy also decrease, since the reservoir 
is then less efficient to counteract decoherence. For pat = 0.2, we get a 1.9 photons state with a 
fidelity of only 54%. Optimizing the average number of atoms per sample is thus important to 
achieve an efficient engineered reservoir. 

Note finally that the phase of the MFSS coherent components is determined by the phase of 
the atomic state superposition when the resonant interaction period begins. Since the atom is 
detuned from the mode during the dispersive interactions, this phase rotates at frequency ±So 
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during the time interval —T/2 < t < —tr/2. The timing of Stark shifts, that determines the 
atom-field interactions, should thus define (T — 1^)/2 with an uncertainty much smaller than I/Sq 
to avoid spurious rotations of this phase. With detunings in the 100 kHz to few MHz range, this 
timing accuracy is easily achieved. 

7 A reservoir for two-mode ESMS 

Our reservoir engineering strategy can be adapted to protect entangled state superpositions of 
two cavity modes, which violate a Bell inequality. Preparation of entangled states of two cavity 
modes, without protection, were considered in [6j [42] . An approximate reservoir for entangling 
large atomic ensembles is proposed and realized in [27]. In ion traps, reservoir engineering has 
recently been used to stabilize a Bell state and a GHZ state of four qubits [26 . 

We present here a scheme in which the two modes belong to the same cavity (two TEM modes 
of orthogonal polarization, whose degeneracy is lifted by an appropriate mirror shape). Extension 
to two separate cavities would require atoms going back and forth between them, a feat not easily 
achieved in the present context of the ENS experiments. 

7.1 Model and target 

We consider two modes a and b of the cavity of respective frequencies uja < oj^. We note b [resp. a] 
the photon annihilation operator for mode b [resp. mode a] and N5 = b'b [resp: N^ = a'^'a] 
the associated photon number operators. A separable joint state of the two modes is written 
li/ja^'ipb)- The atomic qubit (transition frequency cjq ^ Ua^^b) interacts with the modes according 
to the Jaynes-Cummings Hamiltonian, which writes, in a frame rotating at the frequency uJm = 

{ida +CJ6)/2: 

Hjc = A{N,-Na)^^-^{\e){e\-\g){g\) (41) 

+^^(l^)(e|(at+bt)-|e)(^|(a + b)), 

where A = (cj^ — uja)/'^ > and S{t) = co'o(^) —oJm- Here again, S{t) can be adjusted by controlling 
ujQ through the Stark effect. We assume that the couping rt{s) is the same with both modes, a 
restriction that could be easily relaxed. 

We note U the unitary evolution operator associated to Hjc (the overline here denotes two- 
mode operators), solution of the Schrodinger equation: 

^U(t) = -iUjcit) V{t) with U(to) = I . (42) 

We note Ug the two- mode evolution operator corresponding to the parameter set q^ and (M^ '^ , Me '^) 
the associated Kraus operators. Approximate analytical expressions of U^ for the relevant pa- 
rameter sets are given in Appendix |AJ Operators Z and Y generalizing for each mode the ones 
introduced in the previous Sections are also defined in the Appendix. 

Let us consider first the successive resonant interaction of the atoms, intially prepared in 
\uat) = cos(i^/2)|^) + sin(ii/2)|e), with the modes b and a. The corresponding propagator is 

U^ = Y{0^^)Y{O^J. The associated Kraus map ( M^^, M^^ ) stabilizes a tensor product of two 

coherent states |-a,a), where a = 2u/6r for small enough u and Or- 
Under the action of the Kerr-like Hamiltonian 

Hif = -7K((Na+Nb)2 + 2Na) 

for a time tx = 2^~' ^^ initial state \-a^a) would get transformed, up to a global phase factor, 
into: 

\ca) = {\a,a)-i\-a,-a))/V2. (43) 
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Figure 11: Timing of the composite interaction of the atom with the two cavity modes a and b 
at frequencies uj^ > oJa- Bottom frame, sohd hue: time profile of S (difference between the atomic 
frequency ujq and the mean frequency uJm of the two cavity modes) during cavity crossing by one 
atomic sample. For (5 = 0, +A, —A respectively, ujq coincides with Um = ^^^^"^ , cj^, cj^. The tt 
pulse on the atomic state is represented here as a red dot. Bottom frame, dashed line: coupling 
strength Q{vt) with t = when the atom is at cavity center. Top frame: scheme of the propagators 
corresponding to the successive steps in the composite interaction. 

In the next Section, we show that the action of H^ can be simulated by sandwiching the resonant 
interaction U^ between two dispersive interactions. The corresponding reservoir thus stabilizes 

|Ca)- 



7.2 Composite interaction 



The detuning profile S{t) used to stabilize \ca) is represented on Fig. 11 (bottom part). The atomic 
frequency is first set at uJm {S = 0), between t = —T/2 and t = — t^. The atom interacts non- 
resonantly with both modes, with opposite detunings. We restrict in this Section to the dispersive 
regime. The corresponding evolution operator is Z(^(N5 — N^)) (see Appendix [a|, describing 
opposite phase shifts of the two modes driven by the atom, with a phase shift per photon (j). 

The atom is then successively set at resonance with h and a for a time tr. During these short 
time intervals, we neglect the residual dispersive interaction with the other mode. The second 
dispersive interaction with the two modes is performed by setting again (5 = from tr to T/2. 
With this sequence, the phase shifts produced in the dispersive interactions would add up for the 
terms where the atom undergoes |e)(e| and \g){g\ during the resonant parts. Instead, as in the 
single mode case, these phase shifts must cancel out. We thus apply on the atom at t^ a tt pulse 
on the |e) -^ \g) transition. It is driven by a classical source feeding a pulse of microwave with 
negligible duration through the interval between the cavity mirrors. This pulse does not couple 
into the cavity modes. 

The phases of modes a and h evolve at the frequencies ±A. In order to cancel the build-up of 
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these phases during reservoir operation, we constrain the total time T between successive resonant 
interactions to TA = modulo 27r. This condition is easily achieved with the Stark atomic tuning. 
This leads, within irrelevant rotations, to the propagator (see Appendix [A| for a detailed cal- 
culation): 

Ut^U^-" = Z(^(Nb - Na))Y(0j:,J 

Y(0;,JZ(0(N, - Nfc)) . (44) 

Setting the dispersive interactions to produce a (/) = tt phase shift per photon, we get 



eff 



with tK^K — '^/2- The resulting atomic reservoir thus stabilizes the entangled pointer state |cck). 
Adapted detuning profiles 5{t) yield the same propagators when the interaction strength is not 
the same on both modes [43 . Generalization to entangled states with more than two coherent 
components in each mode is straightforward, using slightly more complex detuning sequences. 
Indeed the latter must then be chosen to have additive instead of opposite dispersive effects on 
the two modes. 

7.3 Numerical simulations 



We numerically solve Eq. ( 42 ) and iterate the corresponding Kraus maps starting from the vacuum 
state with u = 7r/4 and 6r = 7r/2, such that the entangled field modes amplitude is of the order of 
1. Decoherence is modeled as the separate coupling of each field mode with a thermal environment, 
with the same damping time Tc and the same temperature {ut = 0.05). The interaction strength 
Q{s) of the atom with each mode has the same Gaussian profile as in the single-mode case, with 
Qo/27T = 50 kHz. In the computations, the field Hilbert space is truncated to the 10 first Fock 
states for each mode. 



Figure 12 shows (solid line) the evolution of the fidelity {ca\p\ca) of the two-mode cavity state 
p w.r.t. an entangled two-component MFSS \ca) with 0.67 photons on the average, starting from 
the vacuum. The reference state is numerically optimized to maximize its fidelity w.r.t. the 
reservoir stationary state (~ p20o)- We have chosen A = Sl^o, ^c = 650 ms. The atomic velocity 
is i; = 22 m/s and each atomic sample has a probability Pat = 0.3 to contain one atom (we neglect 
here two-atom samples). The engineered reservoir is efficient, since the optimal fidelity is ~ 89%. 
This value is reached after ^ 30 samples, corresponding to only 10 atoms on the average. To 
illustrate the protection of the state, we interrupt the reservoir after 200 atomic samples. As 



shown in Figure 12, the fidelity w.r.t. the target state rapidly decreases. 

The entanglement of the state produced by the reservoir can be proved by a violation of a Bell 
inequality adapted to this two- mode case[5[[6]. The Bell signal is: 

B{7a,lb,la,lb) = 'T\Wha,lb) + W{^a,lt) (45) 

+ W^(7;,76)-W(7a,76)|, 

where VK(7a,76) is the two-mode Wigner function. It is defined as: 

W(7a,76) = ^Tr{D%D>L,,p-DlD'^^^P) , 

where P = e*^(^^+^^) is a joint parity operator and D^^ and D^^ are the displacement operators 
for modes a and b respectively. In a local realistic model, B is always smaller than 2. A value 
larger than 2 for some (To^Tb^Ta^Tb) amplitudes is a proof that p is not separable. 



Figure 13 shows a cut of the two-mode Wigner function of p2oo in the plane 5R(7a) = 5R(75) = 
in which maximum violation of the inequality is expected [F. A numerical optimization of the 
Bell signal in this plane provides the four amplitudes shown as white dots. We have performed 
similar optimizations of B after each atomic sample interaction and plotted the maximum Bell 
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Figure 12: Simulation of the reservoir stabilizing a two- mode entangled state. Solid line: fidelity 
of p, the cavity state starting at vacuum, w.r.t. an ideal optimized entangled state of the two 
modes |cq,), as a function of time in units of the sample interaction time T. The reservoir operates 
up to t/T = 200 and is then switched off. Dashed line: maximum Bell signal J3^^^ as a function 
of time. A B^^^ value above the thin dash-dotted line (^J3^^^ = 2) proves entanglement of p. 



signal B^^^ as a dashed line in figure 12 It reaches ^ 2.1 > 2 which implies that the reservoir 
stabilizes a provably entangled state of the modes. When the reservoir is switched off after 200 
interactions, decoherence causes a rapid decrease of B^^^. 



Figure 14 shows the maximum Bell signal B^^^ of the steady state as a function of Tc, for 
three detuning and atomic velocity values. The Bell inequality is violated for all these settings 
when Tc > 450 ms. The crossing of the different curves illustrates the competition between two 
effects. For small TcS, the Bell signal is larger when A is smaller, since a small A corresponds to a 
relatively large velocity and thus to a smaller total interaction duration T. Thus the reservoir is a 
more efficient protection against decoherence when A is small. For very large Tc^ cavity damping 
becomes negligible w.r.t. the dispersive approximation error introduced in the reservoir action, for 
which large A values are preferred. 

The Tc values required for a violation are certainly difficult to reach, but they are only ^ 3 
times larger than the best damping time reported so far [44]. Stabilizing field states violating a 
Bell inequality may thus be within reach of the next generation of experiments. 



8 Concluding remarks 

We have proposed simple engineered reservoirs stabilizing a wide variety of non-classical field 
states in one and two quantum cavity modes. These reservoirs efficiently counteract the standard 
relaxation of the cavities and offer promising perspectives for studies and applications of mesoscopic 
field state superpositions. 

We have gained a detailed insight into the reservoir mechanisms, and particularly into the way 
it corrects for decoherence-induced quantum jumps of the field. We have performed extensive 
numerical simulations justifying the approximations used in [32] and assessing the robustness of 
the method to experimental imperfections. 

We have discussed here, for the sake of definiteness, the reservoir operation in the context of 
the microwave- CQ ED experiments performed with circular Rydberg atoms and superconducting 
cavities at ENS. We have shown that many quantum states protected by our reservoir could 
realistically be observed in this context. Clearly, the method could be straightforwardly extended 
to other spin/spring systems, in cavity QED and trapped ions contexts. It is particularly appealing 
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Figure 13: Cut in the plane (3?(7a) = ^{jb) = 0) of the two-mode Wigner function VK(7a,76) of 
P200- The fringes and negative values for W are a signature of the "quantumness" of the stabilized 
state. The white dots show the points used to maximize the violation of Bell's inequality. 
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Figure 14: Maximum Bell signal J3^^^ of p2oo as a function of the cavity lifetime Tc for A/27r = 
300 kHz, 'u = 30 m/s (solid blue line); A/27r = 400 kHz, 'u = 22 m/s (dashed-dotted green line); 
A/27r = 500 kHz, v = 18 m/s (dashed red line). 

for the thriving field of circuit-QED [34]. Resettable superconducting qubits [33 interacting with 
one or two cavity modes could be used to implement our proposal. With two separate cavities 
interacting with qubits, it would become possible to stabilize a non-local entangled mesoscopic 
superposition and to study the fascinating interplay between decoherence and non-locality. 

A Propagators 

This appendix details the computation of the propagators associated to the atom-cavity interaction 
in the various settings used in the main text. 
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A.l Single- mode case 

For a resonant interaction [S{t) = 0], Eq. Q writes: 

|U(t) = ^(|5)(e|at-|e)(5|a)U(t), 

with s = vt if we set the time origin such that the atom crosses the cavity axis at t = 0. For each 
Bloch sphere Bn associated to the invariant space spanned by (|^, n + 1), |e, n)), this interaction 
induces a Rabi rotation at an angular rate ^/nQ{s) around the Y axis. We therefore define the 
unitary operator Y(/n) [Eq. ([9|] performing a rotation around Y by an angle /(n), where f{n) 
is an arbitrary function of n. The resonant interaction propagator is thus given by Eq. ([l2|. 

For the interaction between ti and ^2 with a constant nonzero detuning S{t) = S y^ [45 , 
the Gaussian variation of Q{vt) precludes an exact integration of Eq. Q. However, assuming 
that Q{vt) varies slowly enough, the coupled atom-field system evolves adiabatically. An initial 
eigenstate of Hjc(ti) (a "dressed state") then remains, for any time t, close to an eigenstate of 
Hjc(t) [4]. This adiabatic approximation is valid provided: 



2^ rse- 



wQ,oyyn-\-l 



2 



(njm) +^"'^'' V.e(^,^), (46) 



for all n in the relevant photon number range. 

The dressed states (| — ,n)^, |+,n)^) that diagonalize Hjc(t) for each n = 1,2, ... satisfy 



and write explicitly 

\-,n)t = cos(ei*V2) |5,n+ 1) + isin(ei*V2) |e,n) 

\+,n)t = isin(^i*V2) |5,n+ 1) +cos(Ci*V2) \e,n) , (47) 

where we define ^n by 

tan^i^) = ^(^t/_2)V^ ^ith^(^)G(^,f). (48) 

The propagator JJq corresponding to the parameter set q = {ti^t2^v^S) is thus 
U, = ^|-,n),,{-,n|,, et<+. 

n 

+ |+,n)i,(+,n|t, e^«+i , (49) 
where the accumulated phase 0^ is given by: 



<Pl = S ^Jl + n{n{vt)/5)^dt. (50) 

Jti 

The restriction of JJq on the Bloch sphere Bn can then be written as: 

= {\-iri)t^{g,n^l\ + \^,n)t^{e,n\) 

x(|^,n + l)(^,n + 1| 6^"^-+! + |e,n)(e,n| e^^^^+i) 

x(|-,n),,(^,n + l| + |+,n),Xe,n|)t . (51) 
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The transformation (| — , n)t(^, n + 1| + |+,n)^(e, n|) is a rotation around the X-axis of Bn 
by an angle — ^^:|.i. The transformation (|^,n + l){g^n + 1| e^^^^+i^^ + |e,n)(e,n| e~'^^^+^^'^) is a 
rotation around the Z-axis of B^ by an angle ^^^i- We thus introduce in Eqs. (|8]),(10) the unitary 



operators X(/n) and Z(/n) representing these rotations on each Bloch sphere Bn. Noting that 
X(— /n)^ = X(/n), we can finally write (|49| in the compact form: 

U,=X(-^(^^))Z(</.?,)X(^(^^)). (52) 

At the start and end of the complete composite interaction, the atom-cavity coupling is weak: 
n^{±vT/2) = 1^2/100. We can thus take X(-^^"^/^^) = X(^^/^^) = I in Section [s] since 
^^^/^^ ~ 0. This leads to Eq. ([25|. 



In the large detuning regime studied in Section |4J we can even neglect all the X operators in 
Uq compared to the large dispersive phase shift operator Z{(j)^). 

A. 2 Two- mode case 

In the two- mode case, it is not possible to get an exact expression for the dressed states. We thus 
restrict either to a resonant interaction with one of the modes or to a dispersive interaction with 
both, assuming a large detuning 2A between modes a and b. In the resonant case, we neglect 
the residual dispersive interaction with the other mode. For the non-resonant interaction, we use 



simple first-order dispersive expressions. In both cases, simulations integrating Eq. (42) explicitly 
confirm the validity of our approximations. 

Let us first investigate the resonant case, with S = ±A. A simple adaptation of the single 
mode results leads to: 

forg = (ti,t2,^,A) (53) 

U, = e-^^(N^-^«)(^2-^^)Z(A(ti-t2))Y(^^J 

forg = (ti,t2,^,-A), (54) 

where Y{0^ ), for instance, is the tensor product of Y(6>^ ) acting on the pair atom-mode a, with 
the identity acting on b. We define a generalized two-mode phase rotation by: 

Z(/n.,nJ = \g){g\ e^ ^---^ + |e)(e| e^ /(-a+D^c-.+D , (55) 

where the operator /Na,Nb is diagonal in the joint Fock state basis of the two modes with 



/Na,Nb l^a^^b) = /(^a,^6)|^a,^6)- In Eqs. (53) and (54), Z is used with a constant argument 
/nJn, =±A(t2-ti). 

We consider now the dispersive interaction corresponding here to J = 0. Applying second-order 
perturbation theory in l^o/^, we get for q = (ti, ^2, 'y, 0): 

U, = e-'^(^^-^")(*^-*^)Z(^'^(Nb-N,)), (56) 

with ^^=2^ f*'n^{vt)dt. 

Using Eqs. ([53|,([54|,(56) and the commutation relation ^, we get an approximate evolution 
operator with the sequence defined in Section [t] (with TA = modulo 27r): 

Ur^uf = V^Z{-A{T/2 + tr)) 

Z(0(Nfc - N,)) Y(^j;,J Y(^j;jJ (57) 

Z(0(Na - Nb)) Z(-A(T/2 - U)) , 

with 



^r\nHvt)dt. (58) 



-T/2 
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The first line in Eq. (57) has no effect on the Kraus map since it is a rotation on the atom only 
after it has interacted with the modes. The operator Z(— A(T/2 — t^)) can simply be compensated 
by properly setting the phase of the Ramsey pulse, preparing now each atom in Z(A(T/2— tr))|iiat)- 



These considerations lead to the effective propagator given in Eq. (44). 



B Equilibrium of reservoir with damping 

If p'^ of the form ([SO]) is a stationary solution of (ISSl) then we have: 



r~ m(^) (/5 - ^z) ((at - z)\z){z\ + |z){z|(a - z)) dz 

J — a'^ 

Kc{p{-z) - /j.{z))z'^\z){z\ dz = , 

(using a.\z) = z\z)^ e*^^a|z) = z\-z) and their Hermitian conjugates). For any real ^, multiplying 
on the left by coherent state (^| and on the right by |^) yields 

2/i(z) {f3 - ^z) (^ - z)e-^^-'^ dz^ KcZ^i^i-z) - /i{z))e-^^-'^ dz = , 

since |(^|2;)p = e~^^~^^ , ^ and z being real. An integration by parts of the first integral yields 
/i(z) (^-^z) e-(^-") - / M| U{z) {p-^z))\ e-(^-") dz 

+ / f^cz'^{/i{-z) - /i(z))e-^^-^^' dz = 0. 



Since this holds for any ^ real, the only possibility is 

Kez\fi{-z) - n{z)) - £ (m(z) (/? - ^z) ) = 

for z G {—o^'oo^^oo) wi^^ ^^^ boundary conditions lim^^c^c^ m(^)(^ ~ <^^) = and /i(— a^) = . 

To solve this differential equation for z G [— Q^^, <^^], we decompose fi in its even and odd parts: 
these parts satisfy two first-order coupled differential equations that can be integrated directly to 



give formula (40) for iii{z 
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